Phase diagram of the dilute magnet LiHo I Yi_ a; F4 
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We study the effective iong-range Ising dipoie modei with a locai exchange interaction appropriate 
for the diiute magnetic compound LiHo I Yi_ I F4. Our calculations yield a value of 0.12 K for 
the nearest neighbor exchange interaction. Using a Monte Carlo method we calculate the phase 
boundary T c (x) between the ferromagnetic and paramagnetic phases. We demonstrate that the 
experimentally observed linear decrease in T c with dilution is not the simple mean-field result, but 
a combination of the effects of fluctuations, the exchange interaction and the hyperfine coupling. 
Furthermore, we find a critical dilution x c — 0.21(2), below which there is no ordering. In agreement 
with recent Monte Carlo simulations on a similar model, we find no evidence of the experimentally 
observed freezing of the glassy state in our calculation. We apply the theory of Stephen and Aharony 
to LiHoa; Yi_ :E F4 and find that the theory does predict a finite-temperature freezing of the spin glass. 
Reasons for the discrepancies are discussed. 

PACS numbers: 75.10.Hk,75.50.Lk,75.40.Mg 
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The rare-earth compound LiHo a; Yi_ x F4 has been 
widely used as a model magnet displaying a wide range 
of phenomena. At T c =1.53 K the predominant long- 
range dipolar interaction causes a second order classical 
phase transition to a ferromagnetic state By applying 
a transverse magnetic field the order can be destroyed in 
a T=0 quantum phase transition at about 4.9 T[2|]. Po- 
sitional disorder can be introduced by substituting the 
magnetic Ho 3+ ions with non-magnetic Y 3+ ions. The 
disorder has been shown to cause a transition to glassy 
behavior at high dilution 

A main attraction of LiHo^Yi_ x F4 is that the micro- 
scopic model is well-known0,|j]. The ground state of the 
Ho 3+ ion in the crystal field is an Ising doublet, with the 
first excited state 11 K above the ground state. At the 
temperature range we consider here (T < 1.5 K) LiHoF4 
should be a very good realization of a dipolar Ising model 
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where J is the dipolar coupling constant, J ex the nearest- 
neighbor exchange constant, the interspin distance 
and Zij the interspin distance along the Ising axis. The 
summation is done over all Ho 3+ ions, which form a 
tetragonal Bravais lattice with four ions per unit cell. 
When diluted, a fraction x of the sites are occupied by 
non-magnetic Yttrium and not included in the above 
sum. The size of the unit cell is (1,1,2.077) in units 
of a = 5. 175 A If we express the interspin distance 
in units of a, then the dipolar coupling constant J = 
(gfj, B /2) 2 /a 3 = Q.214K"0. The exchange coupling J ox 
has been experimentally determined to about half of the 
nearest-neighbor dipolar coupling[5]. In our calculation 
we have neglected the next nearest neighbor exchange in- 
teraction, which was found to be about 5% of the nearest- 
neighbor dipolar coupling^]. In addition, we have left 
out the hyperfine coupling between the nuclear and elec- 



tronic spins as well as the random fields generated by the 
breaking of crystal symmetries due to the dilution. The 
effects of these terms on our results will be discussed. 

A goal of the extensive experimental studies Q of the 
dilute magnet LiHo x Yi_ x F4 is to establish the material 
as a spin glass prototype with canonical glass properties, 
and with a well understood microscopic theory. This 
would allow comparison between different analytical ap- 
proaches to spin-glass systems, as well as provide an im- 
portant experimental benchmark. Currently, it is widely 
believed that the above dipolar Ising model captures the 
essential behavior of LiH02.Y1_2.F4 observed in numerous 
experiments, yet a direct calculation of the phase dia- 
gram is lacking. The goal of this study is to fill this void 
and determine the phase diagram for the dilute dipo- 
lar Ising model appropriate for LiHo j; Yi_ ;r by a direct 
non-approximate Monte Carlo calculation. In the process 
we also address the fundamental question of whether a 
disordered classical dipolar ferromagnet supports a long- 
ranged spin-glass phase. 

The experimentally obtained phase diagram in shown 
in Fig. Q] For x > 0.5 the boundary between the para- 
magnetic and ferromagnetic phases can be fitted to a 
straight line passing through the origin, corresponding 
to the mean-field result T c (x)=xT c (l). As the dilution 
is increased the boundary falls below the mean-field re- 
sult and glassy behavior ensues. At one point (x=0.167) 
freezing of the spin glass was observed and at further di- 
lution (x=0.045) the glassy state did not appear to freeze. 
This so-called anti-glass phase shows a behavior distinct 
from traditional spin glasses and has been the subject of 
numerous investigations 0, 0, 0. 

We are aware of two earlier theoretical investigations 
of randomly parked dipoles. The conclusion of the first 
study H, considering bond-diluted dipoles, was that, de- 
pending on the lattice structure, spin-glass ordering may 
be favored over ferromagnetic ordering at low-T. The or- 
dering (spin glass or ferromagnetic) persists for any finite 
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FIG. 1: Experimental phase diagram from Ref . 0. Open cir- 
cles denote glassy behaviour, SG = spin glass. 



dilution x 1 in disagreement with the anti-glass phase. 
The second study[l(| predicts that a site-diluted BCC 
lattice is ferromagnetically ordered above x=0.21 with a 
spin-glass phase below x=0.21. It is also interesting to 
note that a study of the three dimensional RKKY Ising 
spin glass, with an interaction of mixed sign proportional 
to 1/r 3 , finds that this system lies on the boundary be- 
tween a finite temperature and a T c = spin glass [111]. 

Numerical Monte Carlo studies of dipoles on a dilute 
BCC latticepij find a transition to ferromagnetic order- 
ing at x = 0.3 ±0.1, but are unable to determine whether 
there is a low-T spin glass transition. A more recent 
Monte Carlo study of Ising dipoles on a cubic lat- 
tice at dilutions x=0.045, 0.12 and 0.20 fails to find a 
finite-temperature spin-glass transition. Note that the 
dipolar model on a cubic lattice is not a ferromagnet at 
higher temperatures, unlike L1H0F4. In conclusion, the 
most relevant theoretical and numerical studies to date 
disagree with experiments on the existence and extent of 
the glassy low-T part of the phase diagram. This could 
be partially explained by the subtleties of the dipolar in- 
teraction since numerical and theoretical predictions de- 
pend on the lattice structure and boundary conditions 
used HI, 14 1 . Our goal is therefore to tailor our calcula- 
tions to LiHo a; Yi_ 2 ;F4 in order to be able to compare the 
entire phase diagram with experiments. 

We have studied the dipolar Ising model given by 
Eq. ([1]) using a Monte Carlo method. Due to the long- 
range nature and angular dependence of the Hamiltonian 
this is a challenging problem. Luttinger and Tisza[l]| 
demonstrated that lattice sums depended on the sam- 
ple shape, while Griffiths later showed fl5j] that physical 
properties are independent of sample shape due to break- 
up into sample-shape dependent domains. In LiHoF4 
there is clear experimental evidence for long needle- 



shaped domains jlft 17 1. In order to compare calcu- 
lations to experiments the domain structure has to be 
taken into account, and there are, at present, two differ- 
ent approaches 1J| . Previously the domain structure of 
LiHoF4 was taken into account by performing the Monte 
Carlo simulation over a spherical cavity embedded in a 
cylindrical domain The part of the domain external 
to the cavity is treated in mean-field theory and gives 
rise to an effective field acting on the sphere. 

Here we choose the other approach, which is to impose 
periodic boundary conditions and evaluate the effective 
interaction between spin i and j as a sum over all periodic 
images of spin j. It is important that the thermodynamic 
limit reflects the domain shape. For a needle shaped do- 
main, which is relevant for LiHoF4, this means carrying 
out the sum along the Ising axis prior to the sum in 
the radial direction. A significant speed-up in evaluation 
the sums can be achieved using the Ewald summation 
method, which splits the sum into two rapidly converg- 
ing parts, one in Fourier space, and one in real space. 
The advantages with periodic boundary conditions over 
the cavity method are twofold. The cavity method ne- 
glects all fluctuations outside the spherical cavity while 
the periodic images include at least part of the fluctua- 
tions in the domain. The cavity method was also shown 
to lead to non-monotonic system-size dependence in some 
quantities [1], which is not the case for periodic boundary 
conditions. 

Due to the long-range interactions, the time required 
for one Monte Carlo step scales as N 2 , as opposed 
to N for the short-range case. Adding the compu- 
tational expense of performing disorder averages over 
several hundred copies of the system makes the effi- 
ciency of the Monte Carlo method particularly impor- 
tant. We have therefore compared the efficiency of the 
single spin-flip Metropolis method with continuous time 
Monte Carlo[l8j], the SSE cluster algorithm^ and the 
Wang-Landau method [lj|, which gives explicit access to 
the density of states. In agreement with other studies 
we found that the Wang-Landau method converges very 
slowly for large system sizes. The cluster algorithm al- 
lows for inclusion of a transverse field, but in the present 
low-temperature classical simulations it becomes ineffi- 
cient since all spins tend to join a single cluster. The 
continuous time Monte Carlo method also proved less ef- 
ficient than the traditional single-spin flip, which there- 
fore was used throughout this study. 

In order to determine the extent of the ferromagnetic 
phase, the critical temperature T c is determined as a 
function of disorder x. In the Monte Carlo simulation 
this is accomplished by calculating the Binder ratio for 
the magnetization 
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In addition to the thermal average, an average over 
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quenched disorder configurations d is calculated. The 
critical temperature was extracted from the intersection 
of the Binder ratio for different system sizes. We used 
system sizes up to 10 3 unit cells, containing 4000 spins. 
Disorder averages were performed over a few hundred dis- 
order configurations. A typical run consisted of 2 x 10 6 
Monte Carlo steps of which the first 10 6 steps were dis- 
carded. 

In mean-field theory there are two phases, a 
low-temperature ferromagnetic phase and a high- 
temperature paramagnetic phase separated by a phase 
boundary T c (x) = xT c (l). For the present model T c (l) = 
2.41 K in simple mean-field theory 4], significantly higher 
than the experimental value of 1.53 K. The effects of fluc- 
tuations can be included using a Monte Carlo method, 
and a recent study using the cavity method found that 
T c (l) = 2.03 K[4]. In the present study the periodic 
boundary conditions allow for fluctuations in the do- 
main surrounding the Monte Carlo cell, and we find that 
T c (l) = 1.91 K for the clean system. The difference be- 
tween the present and the experimental result can be 
attributed to an anti-ferromagnetic exchange interaction 
which was measured to about half of the nearest neighbor 
dipolar interaction Treating J cx as a free parameter 
we find that a value of J cx = 0.12 K, or about 38 % of 
the nearest neighbor dipolar interaction J^ ip — 0.33 K, 
lowers T c to 1.53 K. 
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FIG. 2: T c as a function of dilution from experiments (circles) 
and Monte Carlo calculations. The dashed lines represent 
mean-field solutions. 

In Fig. [5] We display the T c (x) boundary for Monte 
Carlo and mean-field theory and compare it to the ex- 
perimental data from Ref. y. 

At low and intermediate dilution, {x < 0.5), the three 
experimental data points follow the mean-field solution. 
In the Monte Carlo data the effects of fluctuations are 
visible already around x — 0.7, particularly without ex- 



change. Including the exchange term makes this effect 
less visible and the Monte Carlo data is in quite good 
agreement with experiments down to x=0.5. However, 
the Monte Carlo data do fall increasingly below the ex- 
perimental results as the dilution is increased. One rea- 
son for this small difference is probably the hyperfine 
coupling between the nuclear and electronic spins[2, 4]. 
This term is important in the low-temperature regime 
and omitted in our analysis. The general effect of the 
hyperfine coupling is to increase the order, and its omis- 
sion would explain why T c (x) decreases faster with higher 
dilution for the Monte Carlo data than for the experi- 
mental data. We have therefore demonstrated that the 
experimentally observed linear decrease in T c is not the 
simple mean-field result, but rather a combination of the 
effects of fluctuations, the exchange interaction and the 
hyperfine coupling. 

In agreement with the experimental data our phase 
boundary appears to intersect the x-axis at a finite value 
of the dilution. This is in sharp contrast to theoretical 
studies [9LIT0II that predict a phase boundary extending to 
the origin. Extrapolating our data the phase boundary 
intersects the x-axis at about x c = 0.15(2) (no exchange), 
and at z c =0.21(2) (including exchange). This is close to 
x = 0.167, where experiments observed freezing of a spin 
glass at T c — 0.13 K. In order to find signs of a spin 
glass freezing we have performed independent simulations 
of two replicas (same quenched disorder) simultaneously 
and the Edwards- Anderson overlap, 

? = £** (1 M 9) . (3) 

i 

has been recorded. For a spin glass freezing to occur 
there should be an intersection of the overlap Binder cu- 
mulants, g q , but no intersection of the magnetic Binder 
cumulant, g m . 

We show the results for the overlap cumulant in Fig. [3] 
The data shown is for the case of no exchange interac- 
tion, but we found similar results when including the 
exchange term. For x = 0.18 the curves intersect around 
T = 0.12 K, but the magnetic Binder cumulant also in- 
tersects at this point, and we conclude that the system 
is magnetized. When we increase the dilution the curves 
do not intersect and we conclude that there is no finite 
temperature freezing of the spin glass above T = 0.05 
K. At temperatures lower than T = 0.05 K equilibration 
problems occur and we cannot exclude the possibility of 
freezing. However, the experimentally observed freezing 
for x — 0.17 occurred at T — 0.13 K, and should be 
visible in our data. 

In order to give further credibility to the phase dia- 
gram in Fig. [2] we plot the magnetization squared as a 
function of disorder in Fig. [4] We note that except for 
the two most diluted systems the finite-size effects are 
very small for the system sizes considered (N=4000 and 
2048). In the limit of high dilution the magnetization 
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FIG. 3: Overlap Binder cumulants in the limit of high dilu- 
tion. 



decreases with increasing system size, indicative of the 
lack of magnetic order. 
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FIG. 4: Magnetization squared for x = n/32 with n = 
4, 8, 12, 16, 20, 24, 28 and 32 (left to right) for N=4000 (dashed 
line) and N=2048 (solid line). 

In order to compare our results to theory we have 
applied the mean-field calculation of Stephen and 
AharonyQ to LiHoF4. The transition temperature for 
the competing ferromagnetic and spin-glass order param- 
eters are given by the two equations 

n = i-J2 xt& M J ij/ k bT c ) = o (4) 

3 

r 2 = 1 -^xtanh 2 (J l3 /fc 6 T c ) = 0. (5) 

3 

For high temperatures r 2 > r\ and ferromagnetic order 



persists, while, depending on the lattice sums, r 2 may 
be smaller than n for low temperatures, in which case 
spin-glass ordering occurs. We have evaluated the sums 
for the lattice appropriate for LiHoF4 and found that the 
solution favors spin-glass order for x c < 0.57. 

One reason for the discrepancy between the experimen- 
tal results and our calculations could lie in parts of the 
Hamiltonian that we have neglected. The hyperfine cou- 
pling between nuclear and electronic spins is important in 
the low-temperature re gime and omitted in our analysis. 
However, a recent study[20( concluded that at zero trans- 
verse field the hyperfine coupling would only renormalizc 
the Ising dipolar Hamiltonian and therefore it should not 
affect the phase diagram qualitatively. In particular, it 
should not be a cause of the spin-glass freezing. Another 
effect omitted in our simulation is the generation of ran- 
dom magnetic fields due to the dilution, which breaks 



the crystalline svmmetrvj20L l2ll . l22j . However, the ef- 
fect of this term should be to increase fluctuations and 
lower the critical temperature for both the ferromagnetic 
and the spin-glass phase. It has even been argued that 
off-diagonal dipolar terms destro y th e spin glass transi- 
tion at any finite transverse field [21j. We conclude that 
not only should the omitted terms not cause a spin-glass 
transition, they also have the potential of destroying the 
long-range glass order. 

The analytic studies [3, [T3| yield the mean-field result 
T c (x) ~ x in the limit of high dilution and therefore pre- 
dict long-range spin glass order extending all the way to 
x = 0. This result differs from both the experimental and 
our numerical studies, which both predict a disordered 
system in the limit of extreme dilution. It therefore ap- 
pears that fluctuations not accounted for in the theory 
are strong enough to cause a finite-dilution phase tran- 
sition at zero temperature. It would be of great interest 
to find a theory that could account for the vanishing of 
the order in the extreme dilution limit. 

Numerical difficulties could also explain the difference 
between our results and experiments. Glassy systems 
are notoriously hard to equilibrate. Energy barriers be- 
tween low-lying states cause equilibration problems and 
make it hard to obtain reliable data for large enough sys- 
tem sizes. The nearest-neighbor Ising spin glass has been 
studied numerically for years, and only recently a consen- 
sus seems to have developed concerning the glass transi- 
tion. In our simulations we see definite signs of equilibra- 
tion problems at the lowest temperatures. In particular 
we find that a decrease in (M 2 ) as the temperature is 
lowered is a clear indicator that the simulation does not 
reach equilibrium. However, having repeated many of the 
simulations we believe that the data we show here is reli- 
able. The system sizes we consider (1000-4000 spins) are 
an order of magnitude larger than in the previous study 
considering dipoles on a cubic lattice [13], but we cannot 
entirely rule out that finite-size effects are so strong in the 
high dilution limit that even larger system sizes would be 
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necessary to see the true thermodynamic behavior of the 
model. 

In order to resolve the differences it would also be im- 
portant to have more extensive experimental data. Wc 
are only aware of two measurements 0, [23| of the spin 
glass transition in LiHoo.i67Yo.833F4- In particular it 
would be of great interest to have further data points 
in the region surrounding x — 0.167 to establish the ex- 
tent and shape of the spin glass phase. Further experi- 
mental data combined with more extensive Monte Carlo 
simulation using parallel tempering, or other improved 
equilibration techniques, should be able to resolve the 
present differences. 
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